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INTRODUCTION 


One  of  the  major  tasks  of  geodesy  Is  the  determination  of  the  gravity  field 
over  the  entire  earth.  Prior  to  the  advent  of  satellites,  global  representation 
of  gravity  *as  severely  hampered  by  sparse  or  poorly  distributed  gravimetric 
measurements  In  the  oceans.  With  the  launching  of  SEASAT  In  1978,  however, 
enormous  Improvements  have  been  achieved  In  our  knowledge  of  the  ocean's  gravity 
field. ^  The  radar  altimeter  readings  from  that  satellite  can  be  used  to  deter¬ 
mine  the  shape  of  the  oceans'  surface  td.th  great  accuracy,  leading  to  important 
applications  not  only  In  geodesy,  but  also  In  geophysics  and  oceanography. 

If  the  earth  were  all  water,  right  to  Its  center,  mean  sea  level  would  be  a 
simple  geometrical  figure  closely  approximated  by  an  ellipsoid  of  revolution. 
However,  the  earth  Is  largely  solid  and  supports  density  variations  in  continen¬ 
tal  areas  and  under  oceans.  These  density  variations  cause  mean  sea  level,  con¬ 
sidered  an  equlpotentlal  surface,  to  vary  by  scores  of  meters  from  an  Ideal  or 
reference  ellipsoid.  The  Irregular,  smooth  surface,  called  the  geold,  formed  by 
the  extension  of  mean  sea  level  throughout  the  earth  is  a  fundamental  surface  In 
geodesy.  Separations  of  the  geold  above  or  below  the  reference  are  dependent  on 
the  variation  of  the  earth's  gravity  field.  This  relationship  has  been  a  prin¬ 
cipal  topic  of  Investigation  by  geodesists,  %fho  have  devised  a  means  for  deter¬ 
mining  the  geold  separations  from  surface  gravimetric  measurements.  Now,  SEASAT 
satellite  altimetry  can  also  be  used  to  derive  geold  separations  and,  thus,  the 
gravity  field  over  the  oceans  with  great  accuracy. 

The  short-arc  technique,  developed  at  the  Air  Force  Geophysics  Laboratory 
(AFGL)  at  Hanscom  AFB,  MA,  provides  us  with  data  representing  the  short  wave¬ 
length  part  of  the  geold;  l.e.,  wavelengths  of  a  few  hundred  kilometers. 

Geoldal  features  at  these  wavelengths  are  primarily  attributable  to  density 
variations  at  deep  levels  In  the  earth's  crust.  The  short -arc  adjustments, 
essentially  a  least-squares  process,  segregates  all  the  short  wavelength  detail 
Into  the  solution  residuals  (geold  separations  or  geold  height  components), 
thus  the  residuals  contain  geold  signals  from  such  sources  as  the  upper  crust 
and  sea -bottom  topography. 

The  short  wavelength  residuals  may  therefore  be  used  to  study  Important 
geophysical  problems  (l.e.  the  mechanism  of  sea-floor  spreading,  crustal 


subductton  at  trenches,  and  global  tides).  One  aspect  that  Is  receiving  atten¬ 
tion  is  the  correlation  of  these  residuals  with  ocean  bottom  topographic  fea¬ 
tures,  particularly  sea  mounts,  whose  accurate  size  and  position  are  required 
for  submarine  navigation. 

A  report^  by  Chat-  -u.  Taljf.ni,  et  al..  Illustrates  many  factors  of  the 
current  research  relating  geold  measurements  In  the  ocean  region  to  gravita¬ 
tional,  density  and  topographic  variations-  However,  a  major  limitation  of  this 
and  prior  studies  was  the  quality  and  quantity  of  the  geold  measurements.  The 
reliability  of  the  input  data  (SEASAT  derived)  enhanced  this  opportunity  to 
better  understand  and  study  the  composition  and  structure  of  the  Earth's  crust 
beneath  the  sea  and  also  the  related  gravitational  variations. 


A.  Calculating  Geold  Anomalies  from  Density  Inhomogeneitles 

The  procedure  here  is  to  calculate  geold  anomalies  for  given  structures  in 
the  earth's  crust.  The  procedure  is  as  follows: 

According  to  Brun's  Formula 

Geold  height  (N)  =  ^ 

where  T  Is  the  anomalous  potential  or  difference  between  the  total  potential  and 
the  normal  potential  due  to  the  reference  ellipse  and  Y  is  normal  gravity.  For 
a  point  mass 


where  G  Is  the  Newtonian  Constant,  M  is  the  mass  and  R  the  distance  between  the 
mass  and  an  observation  point.  Then  for  a  body  (collection  of  point  masses).  In 
Cartesian  coordinates. 


N(x' ,y' 


>■; 


f; 


■ _ p  dx  dy  dz _ 

((x-x')2+(y-y')2+(z.z')2)l/2 


where  p  is  density  and  the  primes  Indicate  the  observation  point. 
A  simpler  two-dimensional  formulation  for  calculating  N  is 


N(x  ,2  ) 


“  '  y  +  (z-z')2)  dx  dz 

IP  one  assumes  that  p  is  sectlonally  constant,  then  this  latter  integral  can  be 
reduced  to  a  line  (contour)  Integral  which  can  be  Integrated  directly.  In 
general,  a  geophysical  model  can  be  comprised  of  Individual  sections  (or 
bodies)  of  constant  density.  For  each  section  the  geold  effect  Is  calculated. 
After  summation  over  all  the  sections  the  total  geold  effect  due  to  the  model  Is 
obtained. 


In  this  study,  the  type  of  section  selected  was  a  two-dimensional  rec¬ 
tangle.  The  geold  height  was  then  calculated  by  both  the  line  Integral  formula 
reported  by  Chapman  and  using  Simpson’s  Rule  for  numerical  Integration. 


z 

Figure  1.  Calculation  of  geold  anomaly  over  a 
two-dimensional  rectangle. 

The  shape  of  the  geold  height  computed  by  both  of  the  above  methods  and 
that  presented  by  Chapman  (Figure  4,  page  82)  are  In  complete  agreement. 

B.  A  Predictive  Algorithm  Relating  Ocean  Topography  to  Geold  Anomalies. 

An  algorithm  to  relate  ocean  topography  to  geold  anomalies  was  developed  by 
relating  bathymetry  (ocean  depth)  to  geoidal  anomalies  In  selected  areas:  speci¬ 
fically,  to  regions  containing  a  significant  oceanic  feature  such  as  a  trench,  a 
sea  mount  or  a  mid-ocean  ridge.  These  regions  are  selected  since  they  are  the 
regions  most  apt  to  show  these  significant  geold  aomalies. 
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The  residual  Capes  were  Chen  ediCed  Co  find  specific  Cracks  In  Che  Hawaiian 
Ridge.  Nine  Cracks  were  selecCed.  The  basic  criCerla  for  selecclon  Is  that 
Che  tracks  be  centered  along  a  significant  topographic  feature.  They  were 
approximately  2600  km  in  length  and  perpedicular  to  Chat  ridge.  These  residuals 
are  considered  only  as  a  function  of  distance  along  the  crack. 

A  program  obcalned  by  AFGL  called  SYNBAPS  was  used  to  provide  bathymetry 
measurements  along  the  same  tracks. 

The  paired  sets  of  data,  residuals  and  bathymetry,  were  then  compared  and 
showed  good  correlation.  The  cross -correlation  coefficients  were  calculated  and 
found  to  be  between  .6  and  .8. 


The  means  of  each  were  then  removed.  No  trends  were  removed  from  the  resi¬ 
duals.  This  In  essence  was  compensated  for  by  Che  fitting  Co  a  geold  model  (GEM 
10).  The  Fourier  Transform  for  each  paired  data  set  was  calculated.  The  ratio 
of  the  paired  transforms  represents  their  admittances. 


admittance 


F(Reslduals) 
F( Bathymetry) 


The  re-transform  of  this  ratio  gives  the  filter  that  when  convolved  with  the 
bathymetry  data  reproduces  the  residuals.  To  determine  one  admittance  or 
filter  function  unique  to  that  area,  the  admittances  or  filters  for  the  nine 
paired  sets  were  averaged  together. 


Simple  physics  suggests  that  the  filter  be  symmetric.  This  means  that  the 
admittance  must  be  real.  This  Is  true  If  the  coherence  =  1,  since  for  each 
wavenumber  of  F(R),  its  phase  is  the  same  as  its  counterpart  In  F(B).  For 
Instance,  let  those  components  be 


F(R)  =  R(R)  +  1  I(R) 


F(B)  =  R(B)  +  1  1(B) 
1  =  , 


then 
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F(R)  ^  L(R1_RJB)_+_.iCR)_i(_b)  ^  J  R(ll  I(R_) 

R(B)2  +  I(B)2  R(B)2  +  I(B)2 

r/b) 

Note  that  the  complex  part  =  0,  If  “  I(By  ’  their  phases  are  equal 

The  Imaginary  part  of  the  averaged  admittance  was  set  to  zero.  The  bathy¬ 
metry  transforms  were  multiplied  by  this  averaged  real  admittance  function.  Re 
transforms  were  then  made  to  obtain  the  filter  In  the  spaclal  domain  and  the 
filter  projected  residuals. 

The  coherence  relating  the  measured  residuals  and  the  bathymetry  was 


obtained  as  follows: 


,  1 

c(k)  “  Fi(R)  Ft(B)* 

i=l 

Ro(R)  =  21  Z  Fi(R)* 


1  f 

Bo(k)  21  ^  ®l(®> 


coherence(k)  ■ 

1  represents  the  number  of  passes  in  the  averaging,  the  *  denotes  the  complex 
conjugate,  ana  ic  represents  the  wavenumber,  R,,  and  Bq  are  the  power  spectra  of 
the  residuals  and  the  bathymetry,  and  c  is  their  cross-spectrum. 

We  then  calculated  the  admittance  Z(k)  from  a  model  for  two  dimensional 
Airy  isostasy: 

Z(k)  -  (Pc-Pw)  e-R^^d-e-kt) 

where 


Pc  density  of  the  crust 
P„  density  of  the  water 
d  is  the  averaged  depth 
t  is  the  thickness  of  crust 
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The  coherence  and  the  comparison  between  the  calculated  and  model  admit¬ 
tance  were  not  at  first  encouraging.  The  bulk  of  the  power  is  frequently 
accounted  for  In  4  or  5  of  the  lower  harmonics.  Noise  and  density  anomalies  are 
expected  in  the  residuals.  These  are  most  apt  to  affect  the  weaker  components 
of  that  Fourier  Transform.  Hence,  the  coherence  for  these  components  becomes  a 
bit  random.  The  caIc  .'  no  a;!;  it;  nee,  being  a  ratio,  also  has  this  problem. 
Hence,  it  does  not  demonstrate  the  decrease  in  amplitude  at  higher  wave  numbers 
as  expected  by  the  reciprocal  of  the  wave  number  term  in  the  Airy  admittance 
model . 

C.  Covariance  Filter  Relating  Bathymetry  Measurements  to  Geoid  Anomalies 

Basically  this  requires  determining  the  coefficients  that  then  can  be  used 
as  a  filter  to  transform  one  type  data  set  Into  another  set.  In  our  case  we  are 
considering  that  filter  as  being  centered  and  symmetric.  Using  the  criteria  of 
Least  Squares  to  determine  these  coefficients,  we  are  then  required  to  minimize 

L-m  m 

I  =  y  [Rs-foBs  -  y  (Bs+i  +  Bs.£)]2 

s=m+l  l=\ 

where  R  and  B  are  paired  data  sets,  f^  the  filter  coefficients,  m=(L-l)/2,  L 
being  the  length  of  the  filter.  Taking  the  partial  derivative  of  I  with  respect 
to  each  fjj,  provides  the  nri-1  normal  (linear)  equations  to  determine  the 
coefficients. 

For  example 

3  T  ^  ™ 

Tf  =  I  [Rs-FoBs-y  F)i  (Bs+Jl+Bs-^)]  (-2)  Bg  =  0 
s=m+l  i!.=l 

^  L-m  m 

aY,  "  L.  [Rs'FoBs’  y  F£  (Bg+Jl+Bg-^)]  (-2)  (Bs+l+Bg-i)  =  0 

s»m+l  )l=l 


o 


m  -  ^  ” 


'J*  "Jl 


^  ^  ^  »•  ^  -■  ••  -•  %•  ."  -"  -  «  *•.  .  •^. 


Residual  and  bathymetry  data  sets  of  approximately  100  points  at  26  km 
spacing  along  specific  tracks  were  obtained.  The  filter  coefficients,  usually 
L»9,  were  determined.  The  filter  was  then  applied  to  the  bathymetry  data  to 
find  "the  predicted  residual  data".  Results  were  reasonable. 


D.  Calculation  of  Gravity  Anomalies  from  Geold  Anomalies. 

In  the  processing  of  the  SEASAT  altimetry  measurement  at  AFGh,  the  geold 
was  considered  as  a  summation  of  the  Standard  Reference  Ellipse  and  the  GEM  10 
model  or  as  It  Is  referred  at  AFGL  the  14-degree  spherical  harmonics  model. 

Their  processing  thereby  resulted  In  separations  (residuals  or  anomalies) 
related  to  that  Input  geold. 

The  following  procedures  were  used  to  calcula  “  the  gravitational  component 
(Ag)  for  these  separations.  The  summation  of  these  Ag  and  the  GEM  10  model  gra¬ 
vity  components  then  represent  the  gravity  anomaly.^  This  anomaly  Is  defined  as 
the  difference  between  measured  surface  gravity  and  that  gravity  computed  for 
the  Standard  Reference  Ellipse. 

First,  we  Introduce  two  formulae  based  on  the  assumption  that  the  oceanic 
geold  separations  are  only  a  function  of  one  dimension  regarding  the  ocean  sur¬ 
face.  This  is  followed  by  two  formulae  which  are  two  dimensional  regarding  the 
ocean  surface. 

The  two  unidimensional  formulae  for  calculating  anomalous  gravity  utilize 
the  separations  as  representing  the  geold  height  or  geold  anomaly  along  an  arc 
a~b  and  are  thereby  related  to  a  corresponding  gravity  anomaly.  The  first  of 
these  formulae^  Is 


Ag(J) 


I  N(b)  _ N(a) 

Ti  a(b)-)l(J)  t(J)-£(a) 


+ 


N(J)-N(X) 

i  (a’(J)-JI(X))*2 


_Y 

Re 


N(J) 


(1) 
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where  Is  the  radius  of  Che  earth,  J  the  point  of  observation  between  a  and  b, 
N()  the  geold  separation  and  the  track  distance.  The  integration  is  to  be 
done  numerically  where  at  the  point  X=J  the  Integrand  is  equated  to  zero. 


The  second  unidimensional  formula  is  derived  from  the  Fundamental  Equation 
of  Physical  Geodesy. ^  Wo  have. 


Ag  =  - 


(2) 


where  r  is  the  direction  normal  to  the  reference  ellipse  surface.  Using 

Y  =  GM/r2  or  --  =  -  —  and 
3f  R 

Brun's  formula  T  =  yN,  relating  the  anomalous  potential  and  geold  height,  we 
find 


Ag  =  - 


3T 

3r 


2yN 

R 


(3) 


where  R  is  the  mean  radius  of  the  reference  ellipse. 


Since  we  are  concerned  with  an  equl potential  surface,  the  surface  of  the 
ocean,  it  is  necessary  for  T  to  satisfy  Laplace's  Equation  in  two  dimensions, 
namely, 


where  z  is  the  direction  normal  to  the  ocean  surface  (approximately  the  same  as 
f ) .  Hence,  we  choose  the  form  of 


E (a|Cos(k£x) 


+  b{Sln(k£x))  e 


-kjz 


where  x  is  distance  along  the  surface  and  kj  Is  the  wave  number.  Using  T/y  “  N 
and  the  distributive  property  of  the  above  equations,  the  separations  are 
thereby  fitted  by  Fourier  series  to  determine  the  aj,  bj.  From  Eq.  (3),  we  find 
at  the  surface  z  =  0, 


(4) 


Ag(x) 


Y  S 
i 


*^1 


2 

R 


(aicosCkjx)  +  bisln(kix)) 


We  now  present  two-dlmenslona'<  surface  formulae  for  calculating  anomalous 
gravity.  If  one  considers  a  thin  circular  ring  centered  about  an  observation 
point  then  one  realizes  that  the  same  geoid  height  measurements  at  different 
locations  around  the  ring  would  contribute  the  same  component  to  the  gravity 
calculation.  One  can  then  closely  estimate  the  anomalous  gravity  due  to  that 
ring  by  using  an  averaged  value  of  geoid  height.  A  concentric  sequence  of  such 
rings  represent  the  neighborhood  of  an  observation  point.  With  this  In  mind, 
and  considering  the  solution  to  Laplace's  Equation  in  cylindrical  coordinates. 
It  Is  noted  that  N  Is  of  the  form 


N(r)  =  C  +  E  aj  J^(rat) 


where  r  represents  the  distance  between  the  point  of  observation  and  point  of 
measured  residual,  Is  the  zero  order  Bessel  Function  and  Is  the  "root" 
where  JoCai)  “0.  If  r  Is  normalized  to  the  distance  to  the  furthest  measure¬ 
ment  considered,  C  Is  the  value  of  N  where  r=l  and 


®1  “  —  [  N(r)  Jo(rai)  rdr  -  — — — r 

J?(ai)  Jq  oiJl(ai) 

The  above  integration  was  done  by  subdividing  the  region  of  Integration  and 
using  averaged  values  for  N(r)  Jolraj)  r  In  each  subdivision.  If  need  be 

b  ^  b 

j  Jo(air)  rdr  *  —  Ji(oir) 

•'a  1  'a 

The  AFGL  data  base  of  observed  separations  are  measured  corrections  to  a 
14 -degree  spherical  harmonics  model.  Hence,  we  need  only  consider  higher  har¬ 
monics  or  that  region  contained  by  a  great  circle  angle  of  (360“/15)/2  or  12° 
to  a  specified  location.  The  units  for  this  wavelength  and  wave  number  are 


Setting 


T  •  Y(C+5:a^jQ(k^r)  e 


where  z=0  corresponds  to  the  earth's  surface,  and  ,  we  find  for  a  select 

position  (z=0,  r*0) ,  that  Eq.  (3)  results  in 

Ag  -  -(Yl-aiki)  -  ^  (C  +  Sai)  (5) 


Next  we  studied  the  double  integral  counterpart  of  Eq.  (1),  namely,  the 
Molodenskli  integral^ 


Ago 


jr 

2x 


(6) 


where  o  designates  the  observation  point  and  1  is  the  distance  between  that 
point  and  an  element  of  surface  do.  Again,  we  computed  averaged  values  for  the 
Integrand  as  a  function  of  distance  along  the  surface  (any  direction)  and  since 
the  area  of  thin  circular  concentric  rings  is  easily  represented,  this  double 
integration  formula  was  evaluated  by  considering  the  integration  replacement 

If ■  ”5  *  -('S---?-.’  • 


The  surface  integration  is  replaced  by  an  inner  circle  and  a  summation  of  con¬ 
centric  rings  centered  about  the  observation  point.  Replacing  the  integrand 
(N-No)/t^  by  its  averaged  value  in  each  ring,  %re  find 


y  N-No  ,  - 

1.1  ’'(r?-r?-i) 


where  (  )i  is  the  averaged  >  For  the  integrand  in  the  i th  ring.  Note  the 
innermost  ring's  contribut.  . 

Examples  of  results  obi  calculations  according  to  formulae  of 

this  section  are  represented  in  l.  following  figures. 


Figure  2.  Profiles  of  gravity  anomaly  surface  measurement  and  that 
obtained  by  using  formula  1. 
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Profiles  of  gravity  anomaly  surface  measurement  and  that 
obtained  by  using  formula  1. 


Figure  3.  ProFiles  of  gravity  anomaly  surface  measurements,  prior  calcula¬ 
tions  made  at  AFGL  using  point  mass  technique  and  that  obtained 
by  using  formula  5  (1*2). 

Figure  4.  A  contour  plot  of  anomalies  obtained  by  using  formula  5  (1=2). 
gravity. 

E.  Auxiliary  Requirements  -  Primarily  Presentation 

In  the  course  of  this  research  effort  a  number  of  sundry  tasks  were 
required.  The  most  notable  task  was  to  represent  AFGL's  deduced  geold 
anomalies. 

The  data  base  for  making  the  contour  plots  of  the  AFGL  measured/computed 
ocean  geold  is  made  from  sorting  their  corresponding  high-pass  filtered  SEASAT 
satellite  altimetry  data.  These  later  data  are  in  sequences  of  points  where 
each  point  contains  the  latitude,  longitude,  14 -degree  geold  model  (GM)  and 
geold  anomaly  or  undulation  (GU).  For  each  of  these  points  the  values  of  GM 
and  GU  is  assigned  to  a  location  in  a  2-dimensional  array.  That  location  is 
determined  by  Indices  obtained  by  rounding  off  the  values  of  latitude  and  longi¬ 
tude  to  the  nearest  Integer.  This  array  represents  a  1  x  1  degree  grid. 

Should  multiple  points  be  assigned  the  same  location,  the  averaged  values  (for 
GM  &  GU)  are  calculated  and  retained. 

The  input  filtered  data  are  restricted  to  ±70”  in  latitude  and  to  ocean 
regions.  The  contour  plotting  program  requires  that  each  location  of  the  array 
have  an  assigned  value.  Should  a  location  not  have  such  a  value,  then  an 
averaged  value  of  the  neighboring  locations  of  ±1  is  substituted.  Except  for 
latitude  =  ±70”  there  are  8  such  neighboring  locations.  The  array  Is  scanned 
for  unasslgned  locations,  and  those  with  8  neighboring  assigned  locations  are 
thus  calculated.  This  Is  then  repeated  for  7,  6,  etc.  This  lowering  of  the 
count  of  assigned  neighboring  locations  Is  needed  to  dummy  fill  the  large  land 
area.  Other  considerations  for  selecting  values  of  GM  and  GU  can  be  optioned  to 
selecting  the  value  of  the  closest  measurement  to  a  grid  point  or  performing  a 
two-dimensional  Interpolation.  The  interpolation  Is  apparently  the  best  but 
also  the  more  difficult,  because  for  each  grid  point  a  value  cannot  be  assigned 
until  all  the  measurements  are  assigned.  Considering  the  number  of  grid  points 
and  the  millions  of  measurements  required,  this  entails  a  rather  Immense 
accounting  problem. 
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Profiles  of  gravity  anomaly  surface  measurements,  prior 
calculations  made  at  AFGL  using  point  mass  tecnnique 
and  that  obtained  by  using  formula  5  (1=2). 


f  anomalous  gravity  obtained  by  using  formula 


The  following  figures  were  selected  as  examples  of  supplemental 
studies. 

Figure  5.  An  Icosohedral  plot  of  the  AFGL  deduced  ocean  geold. 

Figure  6.  A  cartesian  plot  of  selected  satellite  tracks. 

Figure  7.  An  Isometric  plot  of  bathymetry  data. 

F.  Conclusion 

Techniques  have  been  developed  for  the  computation  of  gravity  anomalies 
corresponding  to  geold  anomalies  deduced  from  satellite  altimeter  data.  The 
computations  have  been  done  for  gravity  models  which  are  functions  of  one  and 
two  dimensions  along  the  surface  of  the  ocean.  The  results  are  reasonable. 
Further  Improvement  can  be  expected  when  more  geold  anomaly  data  become 
available. 


cosohedral  plot  of  the  AFGL  deduced  ocean  geoid 


VS'.W 


bathymetry  data 
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